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California  Institute  of  Technology 
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Abstract 

Direct  Numerical  Simulations  (DNS)  of  a supercritical  LOX/i:r2  tempo- 
ral three-dimensional  mixing  layer  are  conducted  for  the  purpose  of  explor- 
ing the  features  of  high  pressure  mixing  behavior.  The  conservation  equa- 
tions are  formulated  according  to  fluctuation-dissipation  (FD)  theory  which 
is  not  only  totally  consistent  with  non-equilibrium  thermodynamics,  but  also 
relates  fluxes  and  forces  from  first  principles.  According  to  FD  theory,  com- 
plementing the  low-pressure  typical  transport  properties  (viscosity,  diffusivity 
and  thermal  conductivity) , the  thermal  diffusion  factor  is  an  additional  trans- 
port property  which  may  play  an  increasingly  important  role  with  increasing 
pressure.  The  Peng-Robinson  equation  of  state  with  a correction  for  obtain- 
ing accurate  molar  volumes,  in  conjunction  with  appropriate  mixing  rules,  is 
coupled  to  the  dynamic  conservation  equations  to  obtain  a closed  system.  The 
boundary  conditions  are  periodic  in  the  streamwise  and  spanwise  directions, 
and  of  non-reflecting  outflow  type  in  the  cross-stream  direction.  Following  the 
DNS  protocol,  the  studied  temperature/pressure  regime  is  one  where  both  Kol- 
mogorov and  Batchelor  scales  can  be  resolved  for  pseudo-species  (i.e.  species 
with  transport  properties  modified  to  allow  the  attainment  of  large  enough 
Reynolds  nmnbers).  Correlations  for  the  Schmidt  and  Prandtl  munbers  as 
functions  of  the  thermodynamic  variables,  based  on  exact  fluid  properties,  are 
used  to  ensure  that  correct  relative  transport  processes  are  employed.  To  ob- 
tain rollup  and  pairing,  the  layer  is  perturbed  similarly  to  heptane/nitrogen 
mixing  layers  that  achieved  transition  in  previous  investigations.  Due  to  the 
strong  density  stratification,  the  layer  is  considerably  more  difiicult  to  entrain 
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than  equivalent  initial  Reynolds  number  gaseous,  droplet-laden  or  supercritical 
heptane/nitrogen  layers.  Simulations  conducted  with  initial  Reynolds  numbers 
of  600  and  750,  an  initial  convective  Mach  number  of  0.4,  an  initial  pressure 
of  100  atm,  and  freestream  temperatures  of  400  K in  the  lower  LOX  and  600 
K in  the  upper  H2  stream,  eventually  exliibit  distorted  regions  of  high  density 
gradient  magnitude  similar  to  the  experimentally  obser\^ed  wisps  of  fluid  at 
the  boimdarj'  of  supercritical  jets.  The  temperature  stratification  was  chosen 
such  that  the  computation  can  be  spatially  resolved  at  these  Reynolds  num- 
bers, accounting  for  memory  constraints.  In  these  simulations  the  layer  does 
not  exliibit  transition  to  turbulence,  although  they  are  conducted  with  an  ini- 
tial Reynolds  number  and  perturbation  for  which  transition  was  obtained  for 
a heptane/nitrogen  mixing  layer  having  a smaller  initial  density  stratification. 
The  cause  of  this  occurrence  is  analyzed  using  global  growth  characteristics, 
vorticity  and  vorticity-niagnitude  budgets,  instcuitaneous  visualizations  of  the 
dynamic  and  thermodynamic  variables  and  an  inspection  of  the  origins  and 
spatial  location  of  dissipation.  The  lack  of  transition  is  traced  to  two  com- 
bined effects.  First,  the  relatively  large  spanwise  perturbation  is  responsible 
for  the  early  creation  of  small  turbulent  scales  which  destroy  the  coherence  of 
the  vortices  formed  through  pairing,  resulting  in  a weakened  ultimate  vortex. 
Second,  regions  of  high  density  gradient  magnitude  are  formed  due  both  to  the 
distortion  of  the  initial  density  stratification  boundary  and  to  the  mixing  of 
the  two  species:  in  these  regions,  the  weakened  vortex  cannot  create  the  small 
turbulent  scales  that  are  crucial  to  mixing  transition. 


1 Introduction 

Liquid  rocket  propulsion  relying  on  hydrogen/liquid-oxygen  (LOX)  combustion  is  not 
a mature  technology  in  that  improvements  in  design  are  still  necessary  to  mitigate 
efficiency  and  stability  problems.  To  solve  these  problems  in  a systematic  manner,  it 
is  required  to  understand  the  fundamental  proeesses  occurring  in  liquid  rocket  cham- 
bers. A simplified  description  of  the  sequence  of  events  in  one  of  these  combustion 
chambers  is  as  follows:  LOX  enters  the  chamber  through  one  of  the  many  injection 
ports,  and  irrespective  of  the  exact  injection  configuration,  LOX  disintegrates,  mixes 
with  Ho  in  a highly  turbulent  manner  while  being  ignited,  with  ensuing  combustion 
producing  water  and  other  incomplete  combustion  products.  From  this  description, 
it  is  immediately  elear  that  LOX  disintegration  plays  a crucial  role  in  determining 
the  size  of  the  LOX  parcels  entering  in  contact  with  H2,  and  further  the  efficiency  of 
the  combustion  process. 

LOX  disintegration  is  a process  essentially  different  from  the  much  studied  spray 
atomization  that  involves  the  breakup  of  a liquid  into  a multitude  of  drops.  Liquid 
breakup  relies  on  physical  mechanisms  involving  the  surface  tension,  and  therefore  it 
is  an  appropriate  concept  only  when  a surface  tension  does  indeed  exist.  In  contrast, 
in  liquid  rocket  chambers  the  mean  pressure  is  about  20  MPa,  with  peaks  as  high 
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as  30  MPa,  and  therefore  both  LOX  and  H2  are  in  a supercritical  state  (see  Table 
1).  By  definition  [1],  a substance  is  in  a supercritical  state  when  it  is  at  a pressure, 
p,  or  temperature,  T,  exceeding  its  critical  value  (indicated  here  by  a subscript  c). 
What  truly  characterizes  the  supercritical  state  is  the  impossibility  of  a two  phase 
region.  Indeed,  when  the  reduced  pressure,  pr  = p/pc  > 1 or  the  reduced  temper- 
ature Tr  = T/Tc  > 1,  in  the  (p,T)  plane  there  is  no  longer  the  possibility  of  a two 
phase  (i.e.  gas/liquid)  region,  and  instead  there  is  only  a single-phase  region  [2].  The 
general  term  for  the  substance  is  fluid,  i.e.  neither  a gas  nor  a liquid.  Noteworthy, 
since  the  critical  locus  of  O2/H2  mixtures  may  include  smaller  or  higher  values  of 
(Pc;  Tc)  than  those  of  the  pure  species,  it  is  possible  that  locally  in  space  and  time 
the  mixture  could  be  at  subcritical  conditions,  but  such  a situation  cannot  be  as- 
sumed to  hold  for  any  amount  of  time  owing  to  the  continuous  change  in  the  mixture 
composition.  Moreover,  Harstad  and  Bellan  [3]  concluded  from  their  study  of  the 
evolution  of  an  initial  LOX  fluid  drop  immersed  in  H2  at  pressures  from  6 to  40  MPa, 
that  nowhere  in  the  fluid  drop  or  in  the  fluid  mixture  surrounding  it,  are  critical 
or  subcritical  conditions  reached.  The  indications  from  this  study  as  well  as  recent 
LOX/i?2  experimental  evidence  from  Mayer  et  al.  [4],  [5]  is  that  LOX  disintegration 
assumes  characteristics  different  from  atomization,  and  that  these  characteristics  re- 
flect features  associated  with  supercritical  conditions.  In  contrast  to  the  process  of 
atomization,  past  the  critical  point  of  the  fluid,  disintegration  assumes  the  aspect 
of  what  Chehroudi  et  al.  [6]  call  ‘fingers’,  or  ‘comb-like  structures’  at  transcritical 
conditions,  having  an  increasingly  gaseous  appearance  with  increasing  pressure;  their 
experiments  were  conducted  with  N2/ N2,  N2/{CO-\-N2),  He/N2  and  02/N2-  Related 
to  the  present  study,  Raman  scattering  measurements  of  the  radial  density  in  free 
N2  jets  at  4 MPa  by  Oschwald  and  Schik  [7]  show  sharp  profiles  independent  of  the 
injection  temperature,  indicating  the  occurrence  of  sharp  density  gradients.  These 
regions  of  sharp  density  gradients  are  indeed  one  of  the  distinctive  optical  features 
in  environments  at  supercritical  conditions.  Not  only  have  they  been  experimen- 
tally observed,  but  they  have  also  been  identified  in  simulations  of  heptane/nitrogen 
three-dimensional  mixing  layers  (see  Miller  et  al.  [8],  and  Okong’o  and  Bellan  [9]). 
Analysis  of  an  enlarged  database  of  Miller  et  al.  [8]  by  Okong’o  and  Bellan  [10]  re- 
vealed that  these  regions  of  large  density  gradient  magnitude  play  a crucial  role  in 
delaying  transition  to  turbulence  by  acting  similar  to  material  surfaces  in  that  they 
damp  emerging  turbulent  eddies. 

Because  the  behavior  of  a binary  species  system  depends  on  the  identity  of  the 
species,  it  is  uncertain  if  our  previous  findings  for  heptane/nitrogen  are  immediately 
applicable  to  the  LOX/i?2  system  which  is  here  of  interest.  For  example,  Harstad 
and  Bellan  [11]  found  that  under  supercritical  conditions,  compared  to  the  hep- 
tane/nitrogen combination,  the  LOX/i72  system  displays  an  increased  solubility,  and 
also  much  larger  effective  Lewis  numbers,  Le^ff-  The  increased  solubility  results  from 
the  thermodynamic  mixing  rules  [1],  whereas  the  enhanced  effective  Lewis  numbers 
were  attributed  to  the  considerably  larger  difference  of  the  specific  (i.e.  mass  based) 
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enthalpies  of  the  components  in  the  LOX/i?2  system  compared  to  the  equivalent 
difference  for  heptane/nitrogen. 

The  present  paper  is  devoted  to  the  study  of  LOX/i/2  three-dimensional  (3D) 
mixing  layers  through  Direct  Numerical  Simulations  (DNS)  as  a means  of  obtain- 
ing information  about  its  specific  behavior.  Since  in  DNS  all  scales  of  the  flow  are 
resolved,  these  simulations  appear  ideal  for  developing  information  to  be  used  in  mod- 
eUng  LOX  disintegration,  as  well  tis  turbulent  LOX/i/2  mixing.  In  Section  2 below 
we  briefly  recall  the  conservation  equations  derived  elsewhere  [8],  [9].  Section  3 is  de- 
voted to  describing  the  particular  aspects  of  the  equation  of  state  implemented  in  the 
model  such  as  to  obtain  increased  accuracy  with  respect  to  the  typical  Peng-Robinson 
formalism.  Further,  in  Section  4 we  address  the  choice  of  the  transport  coefficients. 
The  configuration  and  boundary  conditions  are  addressed  in  Section  5,  whereas  in 
Section  6 we  discuss  the  numerics.  Section  7 focuses  on  the  initial  conditions  and 
presents  results  from  a linear,  inviscid  stability  analysis  which  is  used  to  understand 
specific  aspects  of  the  L0X/i^2  mixing  layer  that  are  necessary  for  choosing  initial 
conditions  for  the  3D  simulations.  In  Section  8,  we  present  two  3D  simulations  at 
different  initial  Reynolds  numbers,  Rcq  . The  global  characteristics  of  the  layers  show 
that  despite  the  relatively  large  momentum  thickness  based  Reynolds  number,  Re^, 
neither  of  these  two  layers  reached  transition.  To  understand  the  origin  of  lack  of 
transition  attainment,  we  concentrate  on  detailed  studies  of  the  layer  with  the  larger 
Reo  and  present  the  results  of  this  analysis.  A summary  and  conclusions  appear  in 
Section  9. 


2 Conservation  equations 


The  conserv^^ation  equations  are  briefly  recalled,  and  the  reader  is  referred  for  details 
to  Miller  et  al.  [8],  Okong’o  and  Bellan  [9]  and  Okong’o  et  al.[12].  For  the  binary 
mixture  under  consideration,  the  conservation  equations  are 


dp  d{puj) 
dt  dxj 


= 0, 


(1) 


djpui)  d {p^ijUj  + pfijj)  ^ dTjj 
dt  dxj  dXj  ’ 


djpYo)  djpYoUj)  djoj 

dt  dxj  dxj  ’ 


djpet)  d [{pet  + p)  Uj]  ^ dgiKj  drijUj 

dt,  dxj  dxj  dxj 
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where  is  a Cartesian  coordinate,  t is  time,  p is  the  density,  Ui  is  the  velocity, 
6t  = e + UiUi/2  is  the  total  energy  (i.e.  internal  energy,  e,  plus  kinetic  energy),  p 
is  the  thermodynamic  pressure  (the  temperature  is  T)  and  Yq  is  the  mass  fraction 
of  O2  (the  mass  fraction  of  H2  is  = 1 - Yo)-  Furthermore,  qiK  is  the  Irwing  - 
Kirkwood  (subscript  IK)  form  of  the  heat  flux  vector  (see  Sarman  and  Evans  [13]), 
jo  is  the  heptane  mass  flux  vector  and  is  the  Newtonian  viscous  stress  tensor 


Tij  = Ai 


dui  ^ duj 
dxj  dxi 


2 duk 
3dxk 


■J%j 


(5) 


where  8ij  is  the  Kronecker  delta  function,  and  p is  the  mixture  viscosity  which  is  in 
general  a function  of  the  thermodjmamic  state  variables.  The  mass  flux  and  heat  flux 
are  given  by 


joj  — 


foj  + — Oih)  YqYh 


pD  dT 
~¥  dxj 


(6) 


QlKj  — 


IK 


dT 

dxj 


OilKRuT- 


m 


morriH 


-fo 


'j’ 


(7) 


j'oj  ~ 


dYo 

dxj 


+ 


YqYh  momH 
RyT  m 


(m 

v,h') 

1 dp' 

\mo 

mn) 

dxj_ 

(8) 


with 


1 momH  ( ^ _ \h 
RyT  m \Tno  mn 


(9) 


The  notation  in  eqs.  6 - 9 is  as  follows:  D is  the  binary  diffusion  coefficient;  ao 
is  the  mass  diffusion  factor  which  is  a thermodynamic  quantity;  rrio,  is  the  molar 
mass  of  species  a;  m = moXo  + itihXh  is  the  mixture  molar  mass  where  the  molar 
fraction  = mYa/nia]  v^a  = {dv / dXa}T,p,Xg{p^a)  is  the  partial  molar  volume  and 
h,a  = {dh/dXa)T,p,Xg{p^a)  IS  the  partial  molar  enthalpy;  v = Xhv^h  + Xqv^o,  is  the 
molar  volume  related  to  the  density  by  n = m/p,  h = Xnh^n  + Xoh^o  is  the  molar 
enthalpy;  Ry  is  the  universal  gas  constant  and  is  a thermal  conductivity  defined 
from  the  transport  matrix 


^iK  = ^ + XhXq  otiK  Q-ekRuRD/ fn, 


(10) 


where  limp_>o  A = Xkt  as  discussed  in  [14],  where  the  subscript  KT  refers  to  Kinetic 
Theory.  The  new  transport  coefficients  associated  with  the  Soret  (in  the  molar  fluxes) 
and  the  Dufour  (in  the  heat  flux)  terms  of  the  transport  matrix  are  asK  and  Qik, 
which  are  the  thermal  diffusion  factors  corresponding  to  the  IK  and  the  Bearman- 
Kirkwood  (subscript  BK)  forms  of  the  heat  flux  (see  Sarman  and  Evans  [13]).  These 
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transport  coefficients  are  characteristic  of  each  particular  species  pairs  and  they  obey 
the  relationship 


— ^IK  ~ OLh-.  (11) 

as  shown  by  Harstad  and  Bellan  [14],  Additionally,  limp_o  ^kt  and  linip-^o  ^bk  — 

afCT- 

To  solve  the  system  of  equations  above,  it  must  be  closed  by  specifying  the  equa- 
tion of  state  (EOS)  and  the  transport  properties. 


3 Equation  of  state 


The  pressure  can  be  calculated  from  the  well-known  Peng-Robinson  (PR)  EOS  given 
the  PR  molar  volume,  vpr,  as 


P = 


RuT 

{^PR  ~ bm) 


i'^PR  + ^bmVpR  - 6?„)  ’ 


(12) 


where  a,„,  bm  are  functions  of  T and  Xa  (see  Appendix  A).  Due  to  the  inaccuracy  of 
the  PR  EOS  at  high  pressures  (see  [1]),  vpr  may  differ  significantly  from  the  actual 
molar  volume  v.  Therefore  for  improved  accuracy,  we  use  a modified  PR  EOS  in 
which  both  vpp  and  the  volume  shift 


Vs  = V-  VpR 


(13) 


are  calculated  from  the  PR  EOS  given  p,  T and  Xa- 

Alt  the  thermodynamic  properties  such  as  the  molar  enthalpy,  h = G—T{dG/dT)p^xt 
the  constant-pressure  molar  heat-capacity,  Cp  — {dh/dT)p^x-,  and  the  speed  of  sound, 

Rs  = i/l//JKs,  are  calculated  in  a consistent  manner  from  the  same  EOS  using  the 
Gibbs  energy,  G: 

G{T,p.Xa)^  [ '^  p{v',T,Xa)dv'  +pv  - R,,T  + J^Xa[Gl  + R^T\n{Xa)]  . (14) 


where  the  superscript  0 represents  the  ‘low  pressure’  reference  condition  for  the  inte- 
gration as  generally  used  in  the  departure  function  formalism  described  by  Prausnitz 
et  al.  [1],  The  integral  is  ill  defined  for  a zero  pressure  reference  condition;  hereinafter 
we  choose  = Ibar  such  that  = RuT/p^.  The  volume  shift  vs  is  calculated  from 
(see  Harstad  et  al.  [15])  as 


The  isentropic  compressibility,  Rg  iss  calculated  from 


(15) 


(16) 
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where  the  expansivity  (qd)  and  the  isothermal  compressibility  {kx)  are  given  by 


(dp/dT),,x  -1 

v{dp/dv)T,x  ’ ^ v{dp/dv)T,x 


(17) 


The  mass  diffusion  factor,  «£),  is  calculated  from  the  fugacity  coefficients,  Pa  (which 
are  related  to  the  Gibbs  energy),  through 


ajD  — 1 + Xa 


d\n{pa) 

dXa 


(18) 


and  portrays  departures  from  mixture  ideality  (i.e.  ao  = !)•  Noteworthy,  ao  is 
independent  of  the  species  chosen  in  the  calculation  of  eq.  18. 

These  equations  specify  the  entire  thermodynamics  of  the  binary  mixture. 


4 Transport  coefficients 


DNS  are  calculations  where  both  the  Kolmogorov  and  Batchelor  scales  must  be  re- 
solved. To  ensure  that  this  requirement  is  satisfied,  we  produced  contour  plots  (not 
shown)  of  the  viscosity,  and  of  Sc  and  the  Prandtl  number,  Pr,  based  on  accurate 
species  transport  properties  calculated  as  in  Harstad  and  Bellan  [16].  Based  on  these 
contour  plots  in  the  range  200K  to  800K,  the  transport  properties  were  correlated  as 


M 


0.75 


(ri  + T2)/2^ 


T in  Kelvins, 


(19) 


Sc  = 


paoD 


= (1.334  - O.668F0  - 0.186Kj  - 0.268FJ) 


1 + 


88.6 


1.51 


(20) 


Pj.  = pCp/m 


1.335 

J'O.l  ’ 


(21) 


where  p,R  is  a reference  viscosity  and  the  reference  temperatures  Ti  (upper  H2  stream) 
and  T2  (lower  O2  stream)  correspond  to  the  free  stream  temperatures  for  mixing  layer 
simulations. 

The  value  of  the  reference  viscosity  is  determined  by  the  specified  value  of  Reo 
(see  below),  chosen  so  as  to  enable  the  resolution  of  all  relevant  length  scales. 

The  thermal  diffusion  factor  is  selected  as  in  Harstad  and  Bellan  [16]  with  asK  = 
0.2,  and  ajx  is  calculated  from  eq.  11. 
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5 Configuration  and  boundary  conditions 


The  temporally  developing  mixing  layer  configuration  is  depicted  in  Fig.  1,  which 
shows  the  definition  of  the  streamwise  (xi),  cross-stream  (3:2)  and  spanwise  (3:3) 
coordinates.  The  layer  is  not  symmetric  in  extent  in  the  X2  direction,  having  found 
in  our  simulations  that  the  layer  gi'owth  is  considerably  larger  in  the  hydrogen  side. 
The  freestream  density  {pi  or  P2)  is  calculated  for  each  pure  species  at  its  freestream 
temperature  (Ti  or  T2)  and  at  the  initial  uniform  pressure  (po)-  The  vorticity  thickness 
is  defined  as  6^{t)  = AC/q/  {d  < ui  > /dx2)^^  where  < > is  the  3:1  — X3  planar 

average  velocity  in  the  streamwise  direction,  and  Af/o  = Ui  — U2  is  the  velocity 
difference  across  the  layer.  Miller  et  al.  [8]  explain  the  choice  of  the  velocities  of  the 
two  streams  in  a simulation  initiated  with  four  streamwise  vortices  pairing  twice  to 
produce  an  ultimate  vortex.  The  choice  of  U\  and  U2  for  a real  fluid 


Ui  — 2Mcfidsi 


(22) 


was  made  with  the  intent  of  keeping  the  ultimate  vortex  stationary  in  the  computa- 
tional domain,  although  the  success  of  this  strategy  was  only  partial.  Here  A/c,o  is  the 
convective  Mach  number  and  Z ~ p/{pTR,Jm)  is  the  compression  factor  indicating 
departures  from  the  perfect  gas  behavior  (i.e.  Z = 1).  The  specification  of  Mc,o 
therefore  determines  AUq.  Given  the  initial  streamwise  velocity  profile  ui  based  on 
Ui  and  Uo,  {d  < Ui  > /d^2)max  hence  = ^ai(O)  ‘‘■re  calculated. 

The  specified  value  of  the  initial  flow  Reynolds  number, 

Re„  = (23) 

PR 


is  used  to  calculate  pr. 

The  boundary  conditions  are  periodic  in  the  streamwise  and  spanwise  directions, 
and  of  outflow  type  for  real  gas  as  derived  by  Okong’o  et  al.  [12].  The  outflow  type 
conditions  are  essential  to  maintain  stability  since  the  initial  perturbation  causes  large 
pressure  waves  which  must  be  allowed  out  of  the  domain  with  minimal  reflection. 


6 Numerics 

6.1  General  method 

The  conservation  equations  are  numerically  solved  using  a fourth-order  explicit  Runge- 
Kutta  time  integration  and  a sixth-order  compact  scheme  for  spatial  derivatives  ([17]). 
Time  stability  is  achieved  by  filtering  the  conservative  variables  every  five  time  steps 
in  the  interior,  in  each  spatial  direction  alternately,  using  an  eighth-order  filter.  Since 
high-order  boundary  filters  are  unstable,  no  filtering  is  applied  at  the  non-periodic 
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{x2}  boundaries.  The  computations  were  parallelized  using  three-dimensional  domain 
decomposition  and  message  passing.  The  tridiagonal  solver  for  the  compact  deriva- 
tive scheme  was  efficiently  parallelized  using  the  method  of  [18].  The  simulations 
were  performed  on  an  SGI  Origin  2000  supercomputer,  using  64  processors. 

In  our  solution  protocol,  once  the  pressure,  temperature,  and  mole  fractions  are 
calculated,  the  density  and  energy  are  calculated  from  the  EOS.  To  calculate  the 
pressure  and  temperature  from  the  known  energy  and  mole  fractions,  we  iterate  at 
each  time  step,  as  described  below. 


6.2  Iterative  solution  for  the  pressure  and  temperature 

Using  the  modified  Peng-Robinson  equation  of  state  [15],  once  (p,  T,  Xa)  are  specified, 
one  can  calculate  (p,  e,  Ya).  However,  in  the  DNS  the  dependent  variables  are  (p,  e,  Ya) 
from  which  v and  Xa  are  calculated  as 


1 mYa 

7Tl  ^ 5 

a.  '^oc 

(y 

rua 


m 

p' 


This  means  that  an  iteration  is  necessary  to  obtain  the  values  of  (p,  T)  corresponding 
to  the  DNS  (p,  e).  The  procedure  for  this  iteration  is  to  update  (p,  T)  at  iteration 
level  n,  (p”,T”),  using: 


jm  ^ j.n-1  ^ pU  ^ pU-L  ^ 


n ^7^  — 1 


dT  = 


j { ^P\  j f ^P\  j f ^P  \ jv 

dp  — I I dv  + I — I de  + „ I dXa- 


Since  Xa  are  known,  dXa  = 0.  Also 


where 


de  = e — e”  dv  = v — 


gn-l  ^ g [p-\T^-\  Xa)  , = V (p^-\ T^-\Xa) 


are  computed  from  the  EOS.  The  needed  derivatives  can  be  calculated  from  Cp,  kt 
and  a„  (where  h,e  in  J/kg,  Cp,Cv  in  J/mol  K)  as 


e,Xa 


CJ  CJ  \de 
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(30) 


(31) 


/ e,Xa  I^T  \ CyV  CyJ  \^^/  v,Xa  Crj  Kt 

The  iteration  is  conducted  by  updating  (p,  T)  until  de  and  dv  are  within  a desired 
tolerance,  e.g.  de/e  < 10“®,  dv/v  < 10“®.  The  initial  guess  for  {p.T)  is  from  the 
previous  time  step,  or  the  previous  Runge-Kutta  stage.  For  the  O2///2  regime  under 
consideration,  (p,  T)  converge  in  2 iterations.  It  was  also  noted  that  if  the  initial  guess 
for  (p;T)  is  (p,c5  (Ti  +T2)/2),  4 or  5 iterations  are  required  for  convergence.  Thus, 
there  appears  to  be  no  need  to  store  (p,  T) , therefore  relaxing  some  of  the  computation 
overhead  associated  with  this  iteration.  In  fact,  the  memory  overhead  of  the  iterative 
method  is  the  addition  of  the  four  arrays  necessary  to  store  e,  v.  e"“\  Compared 
with  the  heptane/nitrogen  simulations  where  a PR  EOS  for  computing  the  pressure 
and  an  energy  fit  for  temperature  were  employed  [9],  here  an  additional  50%  CPU 
time  per  time  step  is  used.  However,  for  the  O2IH2  mixture,  the  form  of  the  energy 
fit  for  the  temperature  would  be  considerably  more  complicated  than  that  of  the 
heptane/nitrogen  mixture,  and  the  pressure  can  no  longer  be  computed  directly  from 
the  molar  volume. 

6.3  Resolution 

The  appropriate  resolution  of  all  scales  is  checked  by  visual  inspection  of  the  dilatation 
field,  V u,  which  is  the  most  sensitive  to  numerical  errors.  The  absence  of  small  scale 
fluctuations  in  V • u is  well  known  to  be  a reasonable  indicator  of  good  resolution. 

Generally,  the  flow  field  is  extremely  sensitive  to  having  an  appropriate  resolution, 
and  its  lack  is  manifested  by  the  code  crashing.  Another  diagnostic  of  inadequate 
resolution  is  an  increasing  number  of  iterations  for  the  convergence  of  the  calculation 
involving  the  EOS,  leading  eventually  to  the  code  crashing  as  well. 

7 Initial  conditions 

The  appropriate  initial  conditions  for  simulating  the  evolution  of  mixing  layers  are 
notoriously  difficult  to  choose,  especially  for  density  stratified  situations  (see  a dis- 
cussion in  Miller  et  al.  [8]).  To  this  end,  following  Drazin  and  Reid  [19],  there  are 
two  issues  that  must  be  addressed:  first,  one  must  inquire  about  the  basic  (i.e.  mean) 
flow,  and  then  about  the  appropriate  disturbance.  Both  of  these  issues  were  thor- 
oughly investigated  by  Okong’o  and  Bellan  [21]  for  real  gases,  and  applied  to  the 
heptane/nitrogen  system  under  supercritical  conditions.  The  same  formalism  is  ap- 
plied here  to  explore  the  specific  features  of  the  Oo/ H2  system.  Having  determined 
here  (by  numerically  solving  the  laminar  equations)  that  the  form  of  the  basic  flow  for 
the  O2/H2  system  is  that  of  an  errorfunction-like  profile  (not  shown),  and  thus  that 
it  has  a single  inflection  point,  a two-dimensional  stability  analysis  is  performed  with 
the  error  function  representing  the  mean  flow;  this  is  acceptable  since  according  to 
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Drazin  and  Reid  [19],  for  mean  flows  displa3dng  a single  inflection  point  the  stability 
analysis  is  not  sensitive  to  the  exact  form  of  the  mean  flow. 

The  freestream  velocity  is  specified  as 

ui{oo)  — Ui;  ui{—oo)  = U2  (32) 

and  the  mean  velocity,  temperature  and  mass  fraction  follow  an  error  function  profile, 
as  listed  below.  Since  the  pressure  is  constant,  the  speed  of  sound  and  inverse  of  the 
density  also  assume  approximately  error  function  profiles.  Therefore,  the  profiles  used 
in  the  stability  analysis  are 

- / X - / X [«1  (oo)  - Ml  (-oo)j  r ./  r-X2\  , J 

(2^2)  — (~oo)  H 2 \ ^ 

a.  (12)  = a.  (-00)  + + ij  _ (34) 

= + [erf(v^|2)+il  (35) 

p{X2)  p(-oo)  2 [p{oo)  p(-oo)J  [ V J 

Several  mean  flow  conditions  are  listed  in  Tables  2-4.  To  find  the  stability  charac- 
teristics of  the  layer,  perturbations  are  imposed  through 

Aui  = Ui  (X2)  exp  [ia  (xi  — ct)] , (36) 

Ap  = p (X2)  exp  [ia  (xi  — ct)] , (37) 

Ap  = p {X2)  exp  [ia  (3:1  — ct)] , (38) 

where  a is  real  and  c is  complex,  and  the  hat  denotes  the  perturbation  amplitudes, 
all  of  which  are  functions  of  X2  only.  Within  the  protocol  of  the  stability  analysis,  the 
physical  quantities  are  obtained  by  taking  the  real  part  (subscript  re)  of  the  complex 
quantities. 

Since  the  mean  velocity,  density  and  speed  of  sound  have  the  same  profiles  as 
for  heptane-nitrogen  [21],  the  stability  curves  could  be  expected  to  be  similar  for 
the  same  initial  density  stratification  and  Mc,o-  However,  as  illustrated  in  Fig.  2a 
where  results  are  obtained  for  the  conditions  listed  in  Tables  2 and  3,  discrepancies 
arise  because,  due  to  the  different  fluid  properties,  the  freestream  Mach  numbers 
are  different.  Generally,  the  most  unstable  wavelengths  are  shghtly  longer  for  the 
O2IH2  mixture,  as  can  be  seen  from  Fig.  2a  and  from  a comparison  of  Table  IV 
in  [21]  with  Table  5 herein,  listing  the  most  unstable  wavelengths.  More  important 
and  directly  relevant  to  the  strategy  of  conducting  simulations  for  temperature  ratios 
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\T2  — Ti  I /T2  as  large  as  possible,  even  for  values  of  | T2  — Ti  | /T-i  quite  smaller  than 
unity,  P2/ Pi  here  very  large,  and  certainly  much  larger  than  for  heptane/nitrogen. 

Equi\’alently,  as  listed  in  Table  I of  [21]  and  Table  3 herein,  P2/P1  = 12.88  corresponds 
to  I T2  — Ti  I /T2  = 0.667  for  heptane/nitrogen  and  0.2  for  O2/H2,  while  j T2  — Ti  | 
/T2  = 0.5  returns  P2/P1  = 24.40  for  O2/H2  as  shown  in  Table  4.  For  smaller  freestream 
temperatures  and  | T2  — Ti  j /T2  = 0.5,  the  stratification  is  even  larger,  as  shown 
in  Fig.  2c  for  a range  of  Ti  and  T2  and  at  p = 100  and  400  atm.  This  indicates 
that  as  P2I Pi  becomes  larger,  due  to  computational  constraints  associated  with  an 
increased  number  of  nodes  in  the  direction  of  the  initial  density  gradient  stratification 
(i.e.  3:2)  for  the  same  resolution,  the  O2/ H2  DNS  must  be  restricted  to  smaller  ratios 
\T2~Ti  I /T2  than  the  equivalent  simulations  conducted  for  heptane/nitrogen  [9], 
[21].  Therefore,  all  simulations  conducted  herein  will  be  for  the  mean  flow  properties 
displayed  in  Table  4,  and  only  Reo  will  be  varied. 

Following  the  arguments  of  Okong’o  and  Bellan  [21]  showing  that  3D  eigensolu- 
tions  to  the  stability  problem  are  not  uniquely  defined,  the  simulations  are  started 
with  heuristic  streamwise  and  spanwise  vorticity  perturbations  [20], [22]  superimposed 
on  the  mean  initial  velocity  profile 


‘^i{x2,x:i)  F^d ^^^-h{x2)h{xi)  (39) 

^ 1 

= F2D^—^fl{xi)f2{x2)  (40) 

1 3 


where  Fi  and  Fa  are  the  circulations, 


We  use  F2D  = 0.1,  Aa  — 1,  Ai  = 0.5,  A2  = A3  = 0.35  for  the  streamwise  perturba- 
tions, and  F30  = 0.05,  Ro  = 1 Bi  = 0.025  for  the  spanwise  perturbations.  The 
wavelength  of  the  perturbation  is  Ai  = 7.296,^  q (the  most  unstable  wavelength  for  in- 
compressible flow)  and  A3  = O.6A1,  following  [20].  Since  this  perturbation  wavelength 
is  smaller  than  the  most  unstable  one,  Ai  = 10.356i^,Oi  and  the  difference  between  the 
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two  values  is  significant,  based  on  the  results  of  Okong’o  and  Bellan  [21]  with  simula- 
tions perturbed  at  wavelengths  smaller  than  the  most  unstable  one,  it  is  expected  that 
even  if  a transitional  state  is  here  achieved,  the  structures  will  not  be  as  convoluted 
as  when  the  layer  would  be  perturbed  with  the  most  unstable  wavelength.  The  initial 
vorticity  thickness  is  6.859x10“^  m.  The  grid  is  chosen  for  all  simulations  so  as 
to  accommodate  four  wavelengths  in  the  streamwise  and  spanwise  directions,  and  the 
evolution  of  the  layer  is  meant  to  encompass  rollup  and  two  pairings  of  the  initial 
spanwise  vortices.  For  these  initial  conditions,  the  grid  sizes  and  the  resolutions  are 
displayed  in  Table  6. 

8 Results 

The  results  of  the  simulations  listed  in  Table  6 are  here  discussed.  For  the  first  sim- 
ulation, the  value  of  Reo  is  chosen  to  be  600  to  emulate  a condition  which  led  to 
a transitional  state  for  the  heptane/nitrogen  mixing  layer  studied  by  Okong’o  and 
Bellan  [9];  however,  that  mixing  layer  had  an  initial  density  stratification  of  12.88 
instead  of  the  much  larger  value  of  24.40  employed  in  the  present  study.  The  larger 
Reo  = 750  of  the  second  simulation  represents  an  effort  to  enhance  the  probability  of 
reaching  a transitional  state;  however,  this  more  elevated  Reo  did  not  lead  to  transi- 
tion either.  To  understand  the  physics  associated  with  these  results,  we  first  focus  on 
the  global  characteristics  of  the  layer  with  special  emphasis  on  the  features  indicative 
of  transition.  Further,  we  investigate  the  specific  aspects  of  the  instantaneous  fields 
which  are  inherently  absent  from  the  measures  given  by  the  global  characteristics. 
Finally,  in  order  to  understand  some  of  the  peculiarities  of  the  instantaneous  fields, 
we  conduct  a second  order  statistical  analysis  focussing  on  the  dissipation. 

8.1  Global  growth  characteristics 

One  of  the  essential  characteristics  of  a mixing  layer  is  its  growth.  Although  many 
definitions  of  growth  appear  in  the  literature,  Cortesi  et  al.  [23]  showed  that  several 
such  measures,  including  the  momentum  thickness,  are  qualitatively  similar.  Here, 
we  define  the  momentum  thickness  as 

2 ^^2, max 

""  7a 'oTi  / (^2+  < pu\  >)(^i+  < pux  >)dx2  (44) 

with  6»1  =<  pui  >x2=i2,ma.  ^nd  02  =<  pui  >x2=L2,mi„»  ^hcre  L2,mm  = -^2/8  and 
T9  may  = 2L2/3.  While  the  growth  is  mostly  a consequence  of  entrainment,  the  prod- 
uct thickness  defined  ^ = f f fy  pYpdV  in  mass  units,  where  Yp  = 2 min(Vo,  Yfj), 

is  a direct  consequence  of  molecular  mixing  as  also  explained  by  Cortesi  et  al.  [23]. 
Both  of  these  quantities,  non-dimensionalized,  are  illustrated  versus  the  nondimen- 
sional  time  t*  = tAUo/S^jfi  in  Fig.  3a,  respectively,  for  the  R600  and  R750  simula- 
tions. The  non-dimensional  momentum  thickness,  8^1  ^wja  of  both  layers  is  similar. 


13 


The  R600  layer  exhibits  the  first  pairing  at  t*  = 80,  but  does  not  show  promise  of  a 
second  pairing  by  t*  = 125,  at  which  time  the  simulation  was  terminated  due  to  lack 
of  further  interest.  In  comparison,  the  R750  layer  pairs  first  at  t*  = 80  and  finishes 
the  second  pairing  at  t*  = 150;  this  simulation  was  continued  to  t*  = 190  to  explore 
the  possible  transition  to  turbulent  mixing  after  the  second  paring.  Despite  the  con- 
tinuous growth  of  the  layer,  and  the  attainment  of  a relatively  large  Re^  = 1680, 
the  R750  layer  does  not  reach  transition  for  reasons  discussed  below.  Compared  to 
other  initially  density  stratified  layers  such  as  drop  laden  layers  (c.f.  Miller  and  Bel- 
lan  [24];  stratification  of  1.5),  or  to  supercritical  layers  of  lesser  initial  stratification 
(c.f.  Okong’o  and  Bellan  [9|;  stratification  of  12.88),  the  present  layers  do  not  show 
the  characteristic  plateaux  indicative  of  the  influence  of  the  forcing.  This 

fact  is  attributed  to  the  much  larger  present  stratification,  resulting  in  a resistance 
to  entrainment.  The  nondimensional  product  thickness,  6p/6p,o,  displays  a continu- 
ous growth  indicating  that  despite  the  lack  of  transition,  molecular  mixing  proceeds 
unabated. 

Depicted  in  Fig.  3b  are  rotational  global  features  of  the  layers;  the  non-dimensional 
positive  spanwise  vorticity,  <<  >>  ^.nd  the  non-dimensional  enstro- 

phy,  <<  ojiUJi  » MJqY , where  u =V  x u is  the  vorticity;  here  <<>>  denotes 
volume  averaging.  Since  the  initial  mean  velocity  profile  is  such  that  the  initial  span- 
wise  vorticity  is  negative,  <<  >>  (^j^  o/At/o)  is  an  indicator  of  the  development 

of  small  turbulent  scales.  Complementing  this  information,  <<  cjjo;,  >>  {8^^/ 
is  a manifestation  of  stretching  and  tilting,  the  mechanism  which  is  primarily  respon- 
sible for  the  formation  of  small  scales.  For  both  simulations,  «oj^  » (<5^,o/ACo) 
increases  from  the  null  value  once  the  layer  rollup  is  completed;  however,  a reduced 
augmentation  rate  is  displayed  by  the  R600  simulation  corresponding  to  the  reduced 
layer  growth.  The  peak  in  the  curve  portraying  the  R600  simulation,  and  the  first 
peak  in  the  corresponding  curve  for  the  R750  simulation  occur  at  the  first  vortex 
pairing.  Another  peak,  but  of  smaller  magnitude  is  displayed  by  the  R750  layer  at 
the  time  station  of  the  second  pairing;  the  relative  magnitude  of  these  two  peaks  is 
a first  indication  of  lack  of  mixing  transition.  Supporting  evidence  of  lack  of  mixing 
transition  evolves  from  examining  the  <<  u/iu^i  >>  (^u;,o/AC/o)^  plots.  The  increase 
in  enstrophy  of  the  R600  is  sporadic  and  modest,  with  two  equivalent  magnitude 
peaks  evident  at  rollup  and  first  pairing;  past  the  first  pairing,  the  enstrophy  decays 
beyond  its  value  at  the  initial  condition.  In  contrast  to  the  R600  results,  for  R750 
the  enstrophy  culminates,  with  a substantial  peak,  at  the  first  vortex  pairing;  sub- 
sequently, the  plot  displays  a decay  supporting  the  lack  of  transition.  According  to 
the  detailed  discussion  in  Okong’o  and  Bellan  [21]  analyzing  the  reasons  for  lack  of 
transition  in  temporal  mixing  layers,  in  the  R750  simulation  we  witness  the  early 
formation  of  substantial  small  turbulent  scales  which  destroy  the  coherence  of  the 
vortices,  thus  impeding  entrainment,  pairing  and  further  development  of  turbulent 
scales. 

The  analysis  presented  below  is  to  ascertain  that  this  physical  picture  is  correct 
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and  complete. 


8.2  Vorticity  production 

To  explore  the  global  rotational  state  of  the  R750  layer,  we  examine  the  budget 
of  the  spanwise  vorticity  in  the  (xi  — 0:3)  homogeneous  planes,  and  inspect  both  the 
average  and  the  RMS.  A parallel  examination  is  conducted  for  the  vorticity  magnitude 
squared. 

The  vorticity  equation  for  a compressible  flow  is 

^ = (a;  • V)u-(V-u)o;-V(-)  x Vp  + V x (^V  • f)  (45) 

Dt  p p 

where  D/Dt  is  the  substantial  derivative,  and  the  equivalent  equation  for  the  vorticity 
magnitude  squared  follows 

=2u-{u-  V)u  - 2(V-u)o;2  - 2u  ■ V(-)  x Vp  + 2a;  • V x (^V  • f).  (46) 

Dt  p p 

The  first  term  in  eqs.  45  and  46  represents  the  stretching  and  tilting  contribution, 
the  second  term  describes  the  eflFect  of  dilatation,  the  third  term  is  the  baroclinic 
participation  to  vorticity  production,  and  the  last  term  portrays  the  viscous  con- 
tribution. Depicted  in  Figs.  4a  and  4b  are  the  non-dimensional  average  and  RMS 
of  the  spanwise  vorticity  budget  at  t*  — 80,  corresponding  to  the  end  of  the  first 
pairing.  Most  of  the  spanwise  rotational  activity  of  the  layer,  both  average  and  RMS, 
is  located  in  the  H2  side,  where  the  lighter  fluid  is  situated.  The  average  spanwise 
vorticity  budget  is  dominated  by  large  peaks  of  the  viscous  term,  while  at  some  lo- 
cations the  stretching  and  tilting  term  competes  with  the  viscous  term.  On  the  LOX 
side  of  the  layer,  the  activity  is  dominated  by  the  barocfinic  term,  while  the  positive 
dilatation  contribution  rivals  in  magnitude  the  negative  stretching  and  tilting  term. 
Compared  to  the  average  spanwise  vorticity,  the  RMS  displays  a large  culmination  of 
the  stretching  and  tilting  term  at  the  boundary  between  the  two  species,  indicating 
that  in  this  crucial  region  of  small  scale  formation  there  is  a considerable  activity, 
explaining  the  large  enstrophy  peak  at  this  time  station.  On  the  H2  side  of  the  layer, 
stretching  and  tilting  and  viscous  terms  contribute  similarly  to  the  RMS  and  dom- 
inate the  dilatation  and  baroclinic  term;  on  the  LOX  side  of  the  layer,  the  viscous 
term  dominates,  although  there  is  appreciable  activity  in  all  other  terms.  The  in- 
dication is  that  in  the  LOX  side,  the  formed  turbulent  scales  are  dissipated  by  the 
action  of  viscosity.  The  result  of  this  dissipation  is  clearly  seen  in  Figs.  4c  and  4d, 
showing  the  non-dimensional  average  and  RMS  of  the  vorticity  magnitude  squared 
budget;  the  viscous  term  is  larger  in  magnitude  than  all  other  terms,  and  negative 
while  the  second  term  in  order  of  decreasing  magnitude  is  the  stretching  and  tilting 
term,  which  is  positive  thus  indicating  production  of  vorticity  by  this  mechanism. 
Production  through  all  mechanisms  is  negligible  in  the  LOX  side  of  the  layer,  and  the 


15 


insignificant  amount  produced  is  dissipated  by  the  dominating  viscous  effect.  Finally, 
the  RMS  of  the  vorticity  magnitude  budget  depicts  the  same  ordering  of  terms  as  the 
RJVIS  of  the  spanwise  vorticity  budget. 

To  understand  the  timewise  evolution  of  the  layer,  illustrated  in  Figs.  5a  and 
5b  are  the  non-dimensional  average  and  RMS  of  the  spanwise  vorticity  budget  at 
t*  = 150.  Similar  to  the  t*  = 80  situation,  most  of  the  spanwise  rotational  activity  of 
the  layer,  both  average  and  RMS,  is  located  on  the  H2  side.  Compared  to  the  mag- 
nitude of  the  equivalent  terms  at  t*  =80,  all  terms  are  here  reduced  approximately 
by  a factor  of  2,  indicating  that  vorticity  production  is  abated;  this  finding  is  totally 
consistent  with  the  global  characteristics  presented  above.  The  dominant  contribu- 
tion to  the  average  spanwise  vorticity  is  from  the  stretching  and  tilting  term  which  is 
negative,  although  at  some  cross-stream  locations  situated  well  into  the  H2  side  of  the 
layer  the  viscous  term  rivals  the  stretching/tilting  one  in  magnitude,  and  is  positive. 
Both  the  averaged  dilatation  and  baroclinic  terms  appear  much  smaller.  This  order- 
ing of  the  relative  magnitudes  is  even  more  dramatic  in  the  RMS  budget.  Clearly,  the 
RMS  production  is  primarily  due  to  stretching/tilting  and  viscosity  effects  which  are 
essentially  of  similar  magnitude;  dilatation  and  baroclinic  influences  are  smaller  by 
approximately  a factor  of  4.  In  Figs  5c  and  5d,  a similar  evaluation  is  presented  for 
the  vorticity  magnitude  squared.  Although  generally  the  same  ordering  of  stretch- 
ing/tilting and  viscosity  versus  dilatation  and  baroclinic  terms  holds,  viscous  effects 
dominate  the  stretching/tilting  activity,  and  the  average  viscous  contribution  seems 
to  extend  further  into  the  LOX  side  of  the  layer.  As  expected,  stretching/tilting  is 
responsible  for  increasing  the  magnitude  of  the  vorticity  squared  due  to  production 
of  small  scales,  whereas  the  viscous  term  diminishes  the  magnitude  of  the  vorticity 
through  dissipation.  The  physical  picture  emerging  is  that  production  of  small  scales 
does  not  keep  up  with  dissipation,  therefore  preventing  the  layer  from  reaching  a 
transitional  mixing  state. 

The  analysis  presented  below  is  devoted  to  (i)  corroboration  of  the  above-derived 
conclusions  based  on  the  global  characteristics  and  vorticity  budgets  of  the  layer  at 
important  time  stations,  (ii)  a documentation  of  the  specific  aspects  of  LOX/H2  mix- 
ing layers  in  the  thermodynamic  regime  chosen  herein,  and  (iii)  an  in-depth  inquiry 
into  the  reasons  responsible  for  the  lack  of  transition.  The  R750  database  is  examined 
exclusively,  but  the  conclusions  pertain  to  the  R600  simulation  as  well. 

8.3  Visualizations  of  the  dynamic  and  thermodynamic  vari- 
ables 

Instantaneous  aspects  of  the  flow  may  reveal  information  that  is  unavailable  from 
global  characteristics.  Such  instantaneous  aspects  are  best  illustrated  through  contour 
plots  at  specific  times  and  locations,  chosen  so  as  to  highlight  important  features  of 
the  layer.  Since  two  important  times  are  t*  — 80  and  150,  as  explained  above,  the 
flow  visualizations  will  depict  the  variables  at  one  or  both  of  those  times. 
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Dynamic  variables  As  discussed  above,  one  of  the  most  fundamental  variables 
indicating  the  evolution  of  the  flow  is  the  spanwise  vorticity.  Shown  in  Figs.  6a  and 
6b  are  braid  cross-section  (0:3  = L3/I6  = 0.0075  m)  distributions  of  the  spanwise 
vorticity  at  t*  = 80  and  150,  respectively;  the  between-the-braid  plane  cross-sections 
(x3  = Z/3/2  = 0.06  m)  display  an  equivalent  behavior.  At  t*  = 80,  the  two  vortices 
remaining  after  the  second  pairing  are  clearly  shown,  whereas  at  t*  — 150  only  a 
single  vortex  appears;  however,  some  remnants  of  the  second  pairing  process  are  ob- 
vious. The  level  of  the  positive  spanwise  vorticity,  indicative  of  small  scale  formation 
decreases  from  t*  = 80  to  150,  consistent  with  the  global  peak  in  positive  spanwise 
vorticity  as  a function  of  time;  however,  the  maximum  positive  spanwise  vorticity  in- 
creases in  the  between-the-braid  plane  as  the  layer  evolves  from  the  flrst  to  the  second 
pairing.  Noteworthy  is  the  irregular,  ‘collapsed’  aspect  of  the  single  vortex,  similar 
to  other  such  single  vortices  resulting  from  two  pairings  in  simulations  that  did  not 
achieve  transition  (see  Okong’o  and  Bellan  [21]).  This  is  due  to  the  early  small  scale 
formation  induced  by  the  relatively  large  value  of  F^Di  and  results  in  the  destruction 
of  coherence  of  the  vortex,  impeding  entrainment  which  is  the  crucial  ingredient  to 
growth. 

The  evolution  of  the  streamwise  vorticity  in  the  mid-braid  plane  was  documented 
by  many  investigators  for  gases  at  atmospheric  conditions,  (e.g.  Rogers  and  Moser 
[25]).  Of  particular  interest  for  qualitative  comparison  with  the  present  results  are  the 
plots  of  Cortesi  et  al.  [27]  for  gravitationally  density-stratified  temporal  mixing  layers; 
for  such  layers,  the  influence  of  gravity  is  measured  by  the  value  of  the  Richardson 
number.  Presented  in  Figs.  7a  and  7b  are  the  present  streamwise  mid-braid  plane 
= 25.6)  contours  at  t*  = 80  and  150,  respectively.  These  instantaneous 
plots  may  be  compared  with  those  in  Figs.  9a  and  9b  of  Rogers  and  Moser  [25] 
corresponding  to  time  stations  after  the  first  and  second  pairings,  respectively.  Both 
in  Figs.  7a  and  7b,  distortions  are  observed  in  the  cross-stream  direction  when 
comparing  with  the  equivalent  figures  in  [25]  which  display  a definite  symmetry.  This 
symmetry  of  the  structures  both  in  the  cross-stream  and  the  spanwise  direction  no 
longer  exists  here  due  to  the  initial  layer  density  stratification.  The  overwhelming 
activity  is  in  the  H2  side  of  the  layer,  containing  the  lighter  fluid.  On  the  other  hand, 
comparisons  with  Fig.  16a  of  Cortesi  et  al.  [27]  show  that  their  streamwise  vorticity 
also  lacks  S5mimetry;  however,  those  calculations  are  not  quantitatively  comparable 
with  ours  because  of  the  different  forcing,  among  other  different  aspects.  For  example, 
Cortesi  et  al.  [27]  find  that  as  the  Richardson  number  is  increased,  the  structures 
recover  some  S5mimetry  (their  Fig.  17a). 

Thermodynamic  variables  In  previous  studies  of  supercritical  mixing  layers  [8], 
[9],  [21]  for  heptane/nitrogen  at  60  atm  and  freestream  temperatures  of  Ti  = 1000  K 
and  T2  = 600  K,  the  peculiarities  of  the  layer  were  associated  with  departures  from 
perfect  gas  behavior  ( Z = 1 for  a perfect  gas),  and  departures  from  mixture  ideality 
(for  an  ideal  mixture  ao  = !)•  These  associations  were  issued  from  quantitative  as 
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well  as  visual  correlations.  Contour  plots  of  Z and  «£>  for  the  present  LOX///2  layer 
(freestream  conditions  of  100  atm,  Ti  = 600  K and  T2  = 400  K)  at  different  time 
stations,  both  in  the  braid  and  in  the  between-the-braid  planes  reveal  that  the  fluid 
is  extremely  close  to  a perfect  gas  and  the  mixture  is  nearly  ideal.  For  example,  at 
t*  = 150  the  compression  factor  varies  between  1.009  and  1.035,  whereas  the  mass 
diffusion  factor  varies  between  0.9939  and  0.9996.  Therefore,  none  of  the  features 
discussed  below  can  be  associated  with  specific  real  gas  or  non-ideal  mixture  behavior. 
It  should,  however,  be  noted  that  this  thermodynamic  state  of  the  layer  could  not  be 
foreseen  a priori,  and  it  is  only  a result  of  the  calculation.  (Moreover,  it  is  expected 
that  at  the  same  pressure,  with  decreasing  temperatures  the  mixture  will  increasingly 
exhibit  departures  from  perfect  gas  and  ideal  mixture  conditions.) 

One  of  the  most  distinctive  features  of  supercritical  mixing  layers  examined  so  far 
is  the  existence  of  regions  of  high  density-gradient  magnitude  (called  further  herein 
by  the  acronym  HDGM).  These  regions  have  been  identified  in  both  pre-transitional 
[8],  [21]  and  transitional  [9],  [21]  supercritical  mixing  layers.  Given  the  perfect  gas 
and  near  ideal  conditions  of  this  layer,  in  retrospect  these  distinctive  features  could 
perhaps  be  better  associated  with  the  initial  density  stratification.  However,  movie 
animations  of  the  [ Vp  ] timewise  evolution  show  that  the  origin  of  these  regions  is  not 
only  from  the  distortion  of  the  initial  boundary  between  the  two  fluid  species,  but  also 
from  the  mixing  between  the  two  species;  this  conclusion  holds  for  all  supercritical 
mixing  layers  studied  so  far,  independent  of  the  binary  system  of  species.  Illustrated 
in  Figs.  8a  and  8b  is  | Vp  ] in  the  braid  and  the  between-the-braid  cross-sections 
located  at  X3  = 0.0075  m and  = 0.06  m,  respectively,  at  t*  = 150.  Compared 
to  heptane/nitrogen  mixing  layers  excited  at  the  same  wavelength  and  F3D  [9],  the 
present  regions  of  high  ] Vp  | are  less  convoluted  and  each  such  structure  is  more 
spread-out,  particularly  in  the  between-the-braid  plane.  The  decreased  convolution 
is  the  result  of  both  lack  of  transition  and  the  fact  that  the  excitation  wavelength  is 
here  further  away  from  the  most  unstable  wavelength  found  in  the  stability  analysis 
(see  Figs.  2a  and  2b).  The  fact  that  each  of  the  HDGM  structures  is  more  spread-out 
is  attributed  to  the  increased  solubility  in  the  LOX/i/2  system  with  respect  to  the 
heptane/nitrogen  one.  For  example,  in  studies  of  heptane  drops  in  nitrogen.  Harstad 
and  Bellan  [14]  found  that  the  initial  density  discontinuity  is  maintained,  although  the 
location  of  the  discontinuity  changes  during  the  drop  evolution.  In  contrast,  a LOX 
drop  in  H2  displays  a quicker  relaxation  of  Yq  from  unity  inside  the  drop,  indicating 
important  H2  solubility  effects  (see  Harstad  and  Bellan  [16]). 

Since  the  existence  of  the  HDGM  regions  cannot  be  here  associated  with  real 
gas  effects,  or  entirely  due  to  the  initial  density  stratification,  the  question  arises  as 
to  their  origin.  Clearly,  they  are  the  result  of  mixing,  without  which  they  would 
not  be  formed.  Parcels  of  heavy  LOX  are  entrained  into  the  lighter  H2  and  they 
lose  their  identity  only  after  mixing  at  the  small  scale.  Before  that  time,  the  much 
larger  molar  weight  of  O2  gives  rise  to  a substantial  density  gradient.  Therefore,  the 
HDGM  regions  are  here  attributed  to  the  very  large  molar  weight  ratio  (a  factor  of 
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16)  between  O2  and  H2-  Contrary  to  the  situation  encountered  for  heptane/nitrogen 
where  diflFusion  (which  is  a molecular  scale  process)  was  inhibited  by  the  lower  than 
unity  mass  diffusion  factor,  here  such  diffusion  is  efficient  {a^  ~ 1).  Thus,  we  find 
that  generally,  HDGM  regions  may  occur  under  quite  different  circumstances.  For 
large  molar  weight  ratio  species,  these  regions  may  occur  even  if  the  species  diffuse 
well  (and  we  note  again  that  this  is  a molecular  scale  process).  For  species  having  a 
modest  to  low  molar  weight  ratio,  HDGM  regions  may  still  occur  if  the  species  have 
difficulty  in  diffusing  into  each  other. 

To  investigate  the  visual  correlation  between  the  HDGM  regions  and  the  distri- 
bution of  Yo,  contour  plots  of  the  LOX  mass  fraction  at  t*  = 150  are  displayed  in 
Fig.  9a  and  9b  representing  the  braid  and  the  between-the-braid  planes.  As  in- 
ferred above,  some  parcels  of  LOX  have  broken  off  from  the  lower  stream  and  are 
seen  convected  to  the  upper  stream.  Also,  the  layer  consists  of  substantially  inho- 
mogeneous fluid  composed  of  both  LOX  and  H2-  An  equivalent  physical  picture  of 
non-homogeneities  is  obtained  when  examining  the  Yq  distribution  at  t*  = 150  in 
the  streamwise  mid-braid  plane  located  at  = 25.6,  shown  in  Fig.  9c.  The 

characteristic  ‘mushroom’  shapes  typical  of  the  streamwise  plane  of  3D  simulations 
are  evident;  equivalent  | Vp  | contours  (not  shown)  are  found.  The  disintegration 
of  the  lower  O2  stream  and  the  migration  of  parcels  of  O2  into  the  H2  stream  are 
manifest.  Noteworthy,  Cortesi  et  al.  [23]  also  detected  these  mushroom  structures 
(their  Fig.  9a)  evolving  during  the  simulation  of  gravitationally  density-stratified  3D 
mixing  layers  forced  deterministically.  These  structures  were  less  developed  when  a 
combined  deterministic-random  velocity  field  was  used  during  the  initialization. 

To  quantitatively  assess  the  composition  of  the  HDGM  regions,  conditional  aver- 
ages are  hsted  in  Table  7 representing  the  spatial  distributions  at  t*  = 150;  a similar 
calculation  performed  at  t*  = 80  yielded  similar  results  up  to  the  third  digit.  The 
conditional  averages  are  calculated  over  regions  of  j Vp  j>|  Vp  |c„to//=  A j Vp  [max 
with  K = 0.1, 0.2, 0.3  and  0.4.  As  higher  values  of  j Vp  j are  probed,  the  mass  frac- 
tion of  LOX  increases,  similarly  to  the  findings  with  the  heptane/nitrogen  system 
[9];  however,  for  the  same  cutoff  constant,  the  average  mass  fraction  values  are  here 
higher.  Nevertheless,  because  of  the  large  ratio  of  O2/H2  molar  weights,  the  aver- 
age Yq  seems  closer  to  stoichiometric  than  could  be  inferred  for  the  heptane/nitrogen 
simulation  for  which  directly  equivalent  evaluations  are  not  possible  owing  to  the  lack 
of  O2  in  the  mixing  layer  simulations.  Considering  that  the  stoichiometric  O2  mass 
fraction  is  here  Yo,s  = 32/(32  -t-  4)  = 0.89,  being  calculated  from  the  reaction 

2H2  + O2 2H2O,  (47) 

the  present  values  of  0.906  to  0.956  found  from  conditional  averages  seem  to  indicate 
a globally  favorable  situation  for  combustion  purposes.  In  contrast,  for  heptane/air 
the  stoichiometric  heptane  mass  fraction  calculated  from  the  reaction 

C^Hiq  -j-  IIO2  (3.7  X ll)iV2  — ^ 7CO2  T 8H2O  -f-  (3.7  x ll)iV2  (48) 
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is  = 100/(11  X 32  + 3.7  x 11  x 28)  = 0.067,  which  compares  less  favorably  with 
values  of  0.84  to  0.921  found  from  similar  conditional  averages  with  K = 0.1. 0.2  and 
0.3  using  the  results  of  a heptane/nitrogen  simulation [9].  This  comparison  is  only 
qualitatively  correct  since  the  simulations  in  [9]  were  based  on  the  heptane/nitrogen 
instead  of  the  heptane/air  system. 

To  further  assess  the  correlation  between  | Vp  | and  Yq,  listed  in  Table  8 are  the 
global  coefficients  found  in  the  braid  and  the  between-the-braid  planes  at  t*  = 150. 
The  volume  based  correlation  between  two  variables  is  defined  by 


« xy  » - «x 


y» 


>/(<<  X-  » - « X >>2)(<<  >>  - <<  T’  >>^) 


(49) 


where  X and  T’  are  generic  variables.  The  correlation  in  Table  8 is  moderate  and 
similar  to  that  found  in  a previous  study  for  heptane/nitrogen. 

Because  the  temperature  is  directly  related  to  the  density  through  the  EOS,  the 
expectation  is  that  the  temperature  distribution  will  be  visually  highly  correlated 
with  I Vp  I if  the  pressure  is  approximately  constant.  Indeed,  examination  of  pressure 
contours  in  both  the  braid  and  the  between-the-braid  planes  (not  shown)  reveals  that 
the  variations  from  the  initial  uniform  pressure  are  small,  at  most  8%.  Consistent 
with  the  almost  uniform  pressure,  the  braid  and  the  between-the-braid  temperature 
contours  shown  in  Figs.  10a  and  10b  are  visually  well  correlated  with  | Vp  |.  Hotter 
fluid  from  the  upper  stream  is  transported  to  the  lower  stream,  and  the  HDGM 
regions  generally  contain  fluid  at  higher  temperature  than  their  surroundings.  This 
situation  is  very  beneficial  to  combustion  since  it  has  been  determined  above  that  in 
these  regions  the  composition  is  close  to  stoichiometric. 


8.4  Irreversible  entropy  production  (dissipation) 

Okong’o  and  Bellan  [10]  have  analyzed  the  reasons  for  the  lack  of  transition  in  3D 
supercritical  heptane/nitrogen  mixing  layers,  and  found  that  it  is  due  to  the  HDGM 
regions  which  acted  similar  to  material  surfaces  and  damped  small  turbulent  scales 
formed  through  stretching  and  tilting.  This  conclusion  was  derived  from  an  irre- 
versible entropy  production  analysis  (the  irreversible  entropy  production  is  in  fact 
the  dissipation).  Specifically,  if  Ej  represents  the  reversible  flux  of  entropy  and  g 
denotes  the  rate  of  irreversible  entropy  production,  then 


Sj  = {qiKj  - liojoj/mo  - fiHjHj/rriH)  /T  (50) 

where  hh  and  po  ^•re  the  chemical  potentials  (partial  molar  Gibbs  free  energy), 
jnj  = —joj-.  2Lnd  g is  the  sum  of  viscous,  Fourier  heat  flux  and  molar  flux  contributions 

Q Qvisc  T gtemp  T 9mass  ; (b^) 
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where  according  to  eqs.  6-8,  gmass  contains  the  departure  from  mixture  non-ideality 
(through  joj)-,  old-,  and  the  Soret  term,  oc  asK-  One  of  the  important  issues  deter- 
mined was  the  contribution  of  each  of  the  terms  listed  in  eq.  52  to  g,  from  (xi  — 0:3) 
plane  averages  considering  both  the  average  and  RMS,  as  well  as  an  equivalent  eval- 
uation from  volume  averages.  Analysis  of  the  data  at  three  time  stations  located 
before,  at  and  after  the  culminating  point  of  the  global  positive  spanwise  vorticity 
showed  that  whereas  at  all  times  the  primary  contribution  to  volume  averages  was 
from  gvisc-  At  the  time  following  the  culminating  global  positive  spanwise  vorticity, 
visualizations  showed  that  most  of  the  g,  as  well  as  the  gmsc  and  g-mass  activity  was 
concentrated  at  HDGM  locations.  This  dissipation  mechanism  was  considered  re- 
sponsible for  the  lack  of  transition,  as  formed  small  scales  were  damped  by  the  region 
of  large  | Vp  | which  acted  similar  to  a material  interface.  These  conclusions  were 
consistent  with  those  of  Hannoun  et  al.  [26]  who  experimentally  investigated  the 
turbulence  structure  near  a sharp  density  interface. 

A similar  analysis  was  performed  by  Okong’o  and  Bellan  [9]  for  a transitional 
supercritical  heptane/nitrogen  mixing  layer,  with  somewhat  different  results.  One  of 
the  important  results  was  that  the  visual  correlation  between  g and  HDGM  locations 
no  longer  existed. 

Similarly  to  our  previous  work,  we  use  here  the  inspection  of  the  dissipation  as  a 
diagnostic  determining  the  causes  of  the  lack  of  transition.  To  this  end,  illustrated  in 
Figs.  11a  and  11b  are  the  braid  and  the  between-the-braid  plane  distributions  of  g at 
t*  =150,  respectively.  A visual  comparison  with  the  | Vp  | plots  of  Figs.  8a  and  8b 
leads  to  the  conclusion  that  there  is  a lack  of  correlation  between  these  two  quantities. 
While  there  are  some  regions  of  common  activity,  few  regions  of  very  large  dissipation 
correspond  to  locations  of  HDGM.  The  regions  of  highest  g and  | Vp  | activity  are 
even  more  separated  at  t*  = 80.  The  difference  between  the  present  situation  and 
that  studied  by  Okong’o  and  Bellan  [10]  for  heptane/nitrogen  is  that  here  the  density 
stratification  is  much  larger,  inhibiting  the  formation  of  small  scales  in  the  LOX  side 
of  the  layer.  This  explanation  is  supported  by  the  vorticity  budget  analysis  conducted 
above  where  the  lack  of  activity  of  the  crucial  stretching  and  tilting  term  on  the  LOX 
side  of  the  layer  was  noted,  and  by  the  striking  visual  lack  of  correlation  between  g 
and  HDGM  locations.  In  contrast,  in  [10]  small  scales  were  formed  on  the  heptane 
side,  but  they  were  being  damped  by  the  action  of  viscosity  in  the  HDGM  regions. 

Homogeneous-plane  average  plots  of  the  g^isc,  gmass  and  gtemp  contributions  to 
g displayed  in  Figs.  12a  (average)  and  12b  (RMS)  at  t*  = 150  are  typical  of  the 
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gvisc  — rp  ^Okk^ll)  , gtemp  ~ J>2  Qy. . . 
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situation  at  t*  = SO  as  well.  Contributions  from  Qmass  and  gtemp  effects  are  aboiit  two 
orders  of  magnitude  smaller  than  viscosity  effects  and  this  holds  for  both  the  average 
and  the  RMS.  Consistent  with  the  discussion  on  the  lack  of  correlation  between  g 
and  the  HDGM  regions  which  are  primarily  located  in  the  LOX  side  of  the  layer, 
here  the  viscous  activity  is  considerably  stronger  on  the  Ho  side  of  the  layer,  where 
the  lighter  fluid  is  located;  this  is  true  for  both  the  average  and  the  RMS.  On  the 
other  hand,  the  dissipation  due  to  the  molar  fluxes  is  stronger  on  the  LOX  side  of 
the  layer  where  the  mass  fraction  and  temperature  gradients  are  larger.  Finally,  the 
much  smaller  gtemp  is  the  result  of  the  enhanced  Leg//  at  supercritical  conditions,  as 
shown  by  Harstad  and  Bellan  [11].  The  larger  Lcg//  causes  the  temperature  to  relax 
faster  than  the  mass  fraction  (non-dimensional  gradients  are  smaller),  as  the  larger 
(than  at  low  pressure  conditions)  thermal  conductivity  promotes  heat  transfer. 

Based  on  this  analysis,  the  lack  of  transition  is  due  to  two  combined  effects. 
Small  scales  are  formed  excessively  early  in  the  evolution  of  the  layer,  destroying  the 
coherence  of  the  ultimate  vortex;  this  was  elucidated  by  inspecting  the  global  char- 
acteristics of  the  layer,  by  evaluating  the  vorticity  and  vorticity  magnitude  budgets, 
and  by  examining  visualizations  of  the  spanwise  vorticity.  As  a result  of  the  weak- 
ened vortex,  entrainment  is  reduced  and  the  small  scales  cannot  develop.  This  effect 
is  enhanced  by  the  very  large  stratification  which  further  prevents  the  formation  of 
small  scales.  Evidence  for  this  latter  effect  is  found  from  exploring  the  development 
of  the  dissipation  and  its  main  contributions,  as  well  as  by  scrutinizing  visualizations 
of  the  dissipation  and  the  density  gradient  magnitude. 


9 Summary  and  conclusions 

A Direct  Numerical  Simulation  study  has  been  performed  of  a three-dimensional  tem- 
poral LOX/H2  mixing  layer  in  order  to  explore  aspects  of  LOX  disintegration  under 
shearing  conditions.  The  conservation  equations  were  based  on  fluctuation-dissipation 
theory,  having  an  enlarged  transport  matrix  that  includes  Soret  and  Dufour  effects. 
To  close  the  system  of  equations,  a real  gas  equation  of  state  was  coupled  to  the 
differential  equations.  Transport  properties  were  accurate  as  much  as  possible  in  the 
chosen  {p,  T)  regime  to  ensure  the  resolution  of  the  Batchelor  scales.  This  was  accom- 
plished by  correlating  the  Schmidt  and  Prandtl  numbers  as  functions  of  the  thermo- 
dynamic variables,  consistent  with  contour  plots  of  these  numbers  based  on  accurate 
properties.  The  values  of  these  numbers  determined  the  thermal  conductivity  and 
diffusivity,  while  the  viscosity  was  determined  from  the  prescribed  initial  value  of  the 
Rejmolds  number.  Boundary  conditions  were  of  periodic  type  in  the  streamwise  and 
spanwise  directions,  and  of  outflow  type,  based  on  a real  gas  characteristic  analysis  of 
the  differential  equations,  in  the  cross-stream  direction.  Additional  to  the  Reynolds 
number,  the  initial  conditions  prescribed  the  Mach  number,  the  temperatures  of  the 
two  freestreams,  the  pressure,  and  the  perturbation  of  the  layer. 

A stability  analysis  conducted  for  the  O2/H2  system  showed  that  at  the  same 
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density  stratification,  the  stability  curve  is  similar  to  the  heptane/nitrogen  system 
previously  studied.  However,  at  the  enlarged  stratifications  characteristic  of  thermo- 
dynamic regimes  of  interest  for  the  O2/H2  system,  the  density  stratification  is  much 
larger  and  the  most  unstable  wavelength  is  longer.  It  is  shown  that  at  high  pressure, 
as  the  oxygen  is  at  lower  temperatures,  the  stratification  increases  dramatically. 

Based  on  the  stability  analysis  and  previous  experience  with  heptane/nitrogen 
supercritical  mixing  layers,  the  perturbation  wavelength  was  chosen  to  be  the  most 
unstable  incompressible  one,  and  the  amplitudes  of  the  excitation  were  those  that 
previously  lead  to  transition  in  a similar  heptane/nitrogen  mixing  layer  simulation 
having  albeit  a smaller  initial  density  stratification.  The  domain  size  was  four  times 
the  perturbation  wavelength  to  accommodate  four  vortices  and  two  pairings.  Two 
mixing  layer  conditions  were  simulated  differing  only  by  the  initial  Reynolds  number. 

As  the  simulations  are  very  computationally  intensive,  one  simulation  was  pursued 
only  past  the  first  pairing,  as  it  was  obvious  that  transition  would  not  be  obtained. 
The  second  simulation,  at  a larger  initial  Reynolds  number,  evolved  through  two 
pairings,  but  also  did  not  reach  transition.  To  determine  the  causes  leading  to  the 
lack  of  transition,  a detailed  analysis  of  the  layer  was  conducted. 

Global  characteristics  of  the  layer  showed  a momentum  thickness  and  product 
thickness  continuous  growth,  with  a relatively  large  momentum  thickness  based  Reynolds 
number  reached.  However,  these  aspects  were  not  sufficient  to  induce  transition.  The 
evolution  of  the  global  positive  spanwise  vorticity  and  the  enstrophy  displayed  a large 
peak  following  the  first  pairing,  and  continued  to  decay  afterwards  with  only  a minor 
increase  following  the  second  pairing.  This  information,  interpreted  in  the  context  of 
a previous  study  examining  the  causes  of  lack  of  transition,  was  a first  indication  that 
the  early  formation  of  small  turbulent  scales  destroyed  the  coherence  of  the  vortices 
formed  after  the  first  pairing,  impeding  entrainment  and  the  further  formation  of 
small  scales. 

To  ascertain  that  this  physical  interpretation  is  correct,  the  vorticity  and  vorticity 
magnitude  budgets  were  scrutinized  at  times  following  each  pairing.  Consistent  with 
the  global  growth  characteristics,  very  fittle  vorticity  is  created  on  the  oxygen  side 
of  the  layer  which  contains  the  heavier  fluid.  Most  of  the  vorticity  is  created  on  the 
hydrogen  side  of  the  layer  by  the  action  of  the  stretching  and  tilting  term.  However, 
even  in  those  regions  the  negative  viscous  term  dominates  the  budget  of  the  vorticity 
magnitude  squared,  draining  vorticity  from  the  system. 

Visualizations  of  the  dynamic  and  thermodynamic  variables  revealed  regions  of 
high  density  gradient  magnitude  which  mostly  exist  in  the  lower  LOX  stream.  These 
regions  are  the  result  of  both  the  distortion  of  the  initial  density  stratification  bound- 
ary and  the  mixing  of  the  two  fluids.  Due  to  the  very  large  molar  weight  ratio  between 
oxygen  and  hydrogen,  parcels  of  LOX  detached  from  the  lower  stream  will  maintain 
their  density  identity  while  being  entrained  into  the  upper,  lighter  hydrogen,  thereby 
creating  these  large  density  gradient  magnitude  regions  prior  to  the  complete  mix- 
ing of  the  two  fluids.  Although  regions  of  high  density  gradient  magnitude  were 
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identified  also  during  the  mixing  of  heptane/nitrogen,  in  that  situation  there  were 
quantitatively  correlated  with  locations  of  mixture  non-ideality.  In  contrast,  at  the 
conditions  of  the  present  simulations  the  fluid  behaves  as  a perfect  gas  and  an  ideal 
mixture. 

Inspection  of  the  irreversible  entropy  production,  which  is  the  dissipation,  con- 
firmed that  most  of  the  activity  is  concentrated  in  the  lighter  upper  stream  where 
most  of  the  small  scales  are  formed.  In  the  lower  stream,  where  most  of  the  high 
density  gradient  regions  reside,  stretching  and  tilting  activity  is  negligible  resulting 
in  the  lack  of  small  scales,  explaining  the  inactive  dissipation.  At  all  locations,  the 
viscous  dissipation  dominates  both  the  mass  flux  and  the  heat  flux  dissipation  by  at 
least  two  orders  of  magnitude. 

According  to  this  analysis,  two  reasons  contribute  to  the  lack  of  mixing  transition. 
First,  the  relatively  large  spanwise  perturbation  induces  early  small  scale  formation 
which  destroys  the  coherence  of  the  vortices  formed  through  pairing,  and  impedes 
entrainment  and  the  further  formation  of  small  scales.  The  ultimate  vortex  resulting 
from  the  second  pairing  is  weakened  during  this  process.  Second,  the  regions  of 
high  density  gradient  magnitude  formed  through  the  distortion  of  the  initial  density 
stratification  boundary  and  also  through  mixing  of  the  two  species  contain  very  dense 
fluid  in  which  small  turbulent  scales  cannot  form  owing  to  the  weakened  ultimate 
vortex. 


10  Appendix  A 

Miscellaneous  relationships  relevant  to  the  EOS  are 

i 3 
i 

where  indices  here  do  not  follow  the  Einstein  notation,  and 

n I r'-'oor*  /r  n 

aij  = 0.45/236 aiaj  (54) 

Pc,ij 


ai  =:  I + Ci  — Ci 


Ci  = 0.37464  + 1.54226Q^  - 0.26992f^f 


bi  = 0.077796^^^ 

Pc,i 
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Tc,ij  = {l-kij)  y/Tc,iTcj  with  ku  = 0 


(55) 


, ..  = iO/3+^;V3y 


^c,ij  — 2 i^c,i  + ^c,j) 


Pc,ij 


RuTc,ijZc, 


with  Tc^i,  Zc,i  Vc,i  and  pc,i  being  the  pure  species  critical  values,  is  the  species 
acentric  factor  and  kij  is  an  empirical  mixing  parameter.  The  values  for  hydrogen 
and  oxygen  are  in  Table  1,  and  for  comparison  the  values  for  heptane  and  nitrogen 
are  hsted  as  well. 

Most  data  references  pertain  to  another  mixing  parameter  related  to 
through 


kij  = 1 


y/audjj 


(56) 


Replacing  in  eq.  56  aij  and  T^^ij  from  eqs.  54  and  55,  yields  a relationship  between 
the  parameters  kij  and  k[j 


(1  - hi)  = (1  - 


c,ij  I 


^c,i3 


1/2 


Given  the  lack  of  information  regarding  the  values  of  k^j  or  k’-p  in  the  simulations 
herein  kij  = 0. 
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Species  m (g/mol) 

Tc  (K)  Pc(MPa)  n, 

;(cm'Vmol)  Zc 

0 

Hydrogen  2.0159 

33.0  1.2838 

64.284  0.306 

-0.216 

Nitrogen  28.013 

126.26  3.399 

89.8  0.290 

0.039 

Oxygen  31.9988 

154.58  5.0430 

73.368  0.288 

0.025 

Heptane  100.205 

540.2  2.74 

432.0  0.263 

0.35 

Table  1;  Pure  species  properties. 

Mean  quantity 

ajo  = — oo  (Oxygen) 

X2  = oo(Oxygen) 

ui(m/s) 

-193.732 

193.732 

as(m/s) 

484.329 

484.329 

p(kg/m-‘’) 

63.191 

63.191 

p(atm) 

100 

100 

T(K) 

600 

600 

1 

1 

Table  2 

: Mean  flow  properties  p2/pi=l- 

Mean  quantity 

X2  = —oo  (Oxygen) 

X2  = oo(  Hydrogen) 

«i(m/s) 

-187.287 

666.798 

as(m/s) 

467.045 

1671.193 

pCkg/m"*) 

68.271 

5.298 

p(atm) 

100 

100 

T(K) 

556 

444 

lo 

1 

0 

Table  3:  Mean  flow  properties  po/pi=12.88. 
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Mean  quantity 

X2  = —00  (Oxygen) 

X2  = oo(Hydrogen) 

Mi(m/s) 

-158.004 

770.983 

as(m/s) 

397.517 

1915.376 

p(kg/m^) 

96.764 

3.965 

p(atm) 

100 

100 

T(K) 

400 

600 

Yo 

1 

0 

Table  4:  Mean  flow  properties  p2/pi=24.40. 


Case 

Flow  type 

El 

Pi 

c 

Cj-e 

~u 

II 

1 

Variable  density^  {as  = 10®) 

1.00 

0.860 

0.3830 

0.6598 

7.31 

2 

Variable  density^ 

1.00 

0.797 

0.3151 

0 

7.88 

3 

Variable  density^ 

12.88 

0.670 

0.1756 

-0.0747 

9.38 

4 

Variable  density^ 

24.40 

0.607 

0.1284 

-0.0745 

10.35 

^Same  velocity  profile 

as  for 

Case  4. 

^Velocity  profile  from  Equation  32,  Mf.fi  = 0.4. 

Table  5:  Most  unstable  wavelength,  two-dimensional  analysis. 


Run  Reo  Li  x L2  x L3  N\  y.  N2  y.  N3  max  Re^^^  tmax  Timesteps  CPU(h) 
R600  600  0.2x0.232x0.12  288x336x176  1014  127.45  3730  5472 

R750  750  0.2  x0.2  x0.12  352  x 352  x 208  1680  190.44  6860  13214 

Table  6;  Listing  of  the  simulations  and  associated  resolution.  Li  is  in  meters. 


1 Vp  \cutoff — 
K 1 vp  Imax 

Braid  plane: 

1 Vp  Uax=  2.440  X lO^kg/m'^ 

Between-the-braid  plane: 

1 Vp  |max=  1-955  X lO^kg/m'^ 

K = 0.4 

0.954 

0.956 

K = 0.S 

0.945 

0.951 

K = 0.2 

0.930 

0.949 

K = 0.1 

0.906 

0.923 

Table  7:  Conditional  averages  over  regions  where 

Vp  |>cutoff.  The  calculations 

made  at  t*  = 150  for  R750. 


Braid  plane  Between-the-braid  plane  Global 

Yo,  I Vp|  Correlation  0.34  0.37  0.35 

Table  8:  Correlations  with  ] Vp  | at  = 150  for  R750. 
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b) 


Figure  2:  a)  Comparison  of  the  stability  curves  for  heptane/nitrogen  (at 
60  atm)  and  oxygen/hydrogen  (at  100  atm),  b)  stability  curves  for  oxy- 
gen/hydrogen (at  100  atm)  and  c)  density  ratio  versus  temperature  for  oxy- 
gen/hydrogen. 
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Figure  3:  Time  evolution  of  a)  momentum  and  product  thicknesses  and  b) 
global  spanwise  vorticity  and  enstrophy. 
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Figure  4:  Vorticity  budget  for  R750  at  t*=S0:  Spanwise  vorticity,  a)  plane 
average  and  b)  plane  RMS,  and  vorticity  magnitude,  c)  plane  average  and 
d)  plane  RMS. 


d)  ^^^(0,0 


Figure  4:  (continued)  Vorticity  budget  for  R750  at  t*=80:  Spanwise  vorticity, 
a)  plane  average  and  b)  plane  RMS,  and  vorticity  magnitude,  c)  plane  average 
and  d)  plane  RMS. 
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Figure  5:  Vorticily  budget  for  R750  at  f*=150:  spanwise  vorticity,  a)  plane 
average  and  b)  plane  RMS,  and  vorticity  magnitude,  c)  plane  average  and 
d)  plane  RMS. 
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Figure  5:  (continued)  Vorticity  budget  for  R750  at  ^=150;  spanwise  vortic- 
ity,  a)  plane  average  and  b)  plane  RMS,  and  vorticity  magnitude,  c)  plane 
average  and  d)  plane  RMS. 
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Figure  6:  Nondimensional  spanwise  vorticity  for  R750  in  the  braid  plane 
(x3  = L3/I6)  at  a)  r=80  and  b)  r=150. 
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Figure  7:  Nondimensional  streamwise  vorticity  for  R750  in  the  streamwise 
mid-braid  plane  (rri/^^,o=25.6)  at  a)  r=80  and  b)  f*=150. 
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Figure  9:  (continued)  Oxygen  mass  fraction  for  R750  at  r=150,  a)  in  the 
braid  plane,  b)  in  the  between-the-braid  plane  and  c)  in  the  streamwise  mid- 
braid plane. 
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Figure  11:  Dissipation  (in  J/m^K)  for  R750  at  f*=150,  a)  in  the  braid  plane 
and  b)  in  the  between-the-braid  plane. 
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Figure  12:  Contributions  to  the  dissipation  (in  J/m^K)  for  R750  at  r=150, 
a)  plane  average  and  b)  plane  RMS. 


